home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Collection of Tools & Utilities
/
Collection of Tools and Utilities.iso
/
basic
/
chaosexe.zip
/
XPHASE3D.TRU
< prev
next >
Wrap
Text File
|
1980-01-01
|
3KB
|
94 lines
!PROGRAM PEND/PHASE3D
LIBRARY "3dlib.trc"
CLEAR
PRINT" ***PENDULUM - 3D PHASE DIAGRAM***"
PRINT
PRINT"THIS PROGRAM SHOWS THE COMPLETE 3D PHASE SPACE OF THE PENDULUM. THE"
PRINT"TIME COORD IS 3/2 OF THE CYCLIC PHI COORD, AND THEREFORE THE BOUNDARY"
PRINT"CONDITIONS ON TIME ARE ALSO CYCLIC."
!
declare def accel
input prompt"Input driving force strength: ":g
input prompt"Input damping :":q
input prompt"Input initial position: ": xint
input prompt"Input initial angular velocity:":v
input Prompt"Input min. and max. time:":tmin,tmax
input prompt"Input 'cyclic' time step as fraction of phi: ":phistep
let w=0.6666667
let eps=1.0e-8
let tstep=0.001
let t=0
let x=xint
Call Parawindow(0,3*pi,-pi,pi,-3,3,s$)
Set mode "hires"
Call Ticks3(1,1,1,s$)
CALL PLOTTEXT3(0,2,5,"PENDULUM - 3D PHASE DIAGRAM",S$)
CALL PLOTTEXT3(0,5,2,"G="&STR$(G)&" Q=2",S$)
let phistep=2*pi*phistep
for i=1 to 1000000
call rk4(x,v,tstep,xnew,vnew,t,w,g,q)! Take a step of size tstep
let tshalf=tstep/2
call rk4(x,v,tshalf,xnh,vnh,t,w,g,q)!Take two half steps
call rk4(xnh,vnh,tshalf,xn,vn,t+tshalf,w,g,q)
let d1=abs(xn-xnew)
let d2=abs(vn-vnew)
let delta=max(d1,d2)
if delta<eps then
if t>tmin then
let tnew=t+tstep
let w1=mod(-w*t,phistep) !Check for time axis step
let w2=mod(w*tnew,phistep)
if w1<w*tstep then
if w2<w*tstep then
let ts=w1/w
call rk4(x,v,ts,xp,vp,t,w,g,q)
if abs(xp)>pi then let xp=xp-2*pi*abs(xp)/xp
let period=2*pi/w
let tyme=mod(t+ts,period)
Call PlotOff3(tyme,xp,vp,s$)
end if
end if
end if
let x=xnew
let v=vnew
let t=t+tstep !Expand step size
let tstep=tstep*.95*(eps/delta)^.25
if abs(x)>pi then !bring theta back into range
let x=x-2*pi*abs(x)/x
end if
else !else reduce step size
let tstep=tstep*.95*(eps/delta)^.2
end if
if t>tmax then let i=1000001
next i
call plottext3(10,0,0,"TIME",s$)
call plottext3(0,4,0,"ANGLE",s$)
call plottext3(0,-1,4,"ANG. VEL.",s$)
end
!
!
sub rk4(x,v,tstep,xnew,vnew,t,w,g,q)
declare def accel
let xk1=tstep*v
let vk1=tstep*accel(x,v,t,w,g,q)
let xk2=tstep*(v+vk1/2)
let vk2=tstep*accel(x+xk1/2,v+vk1/2,t+tstep/2,w,g,q)
let xk3=tstep*(v+vk2/2)
let vk3=tstep*accel(x+xk2/2,v+vk2/2,t+tstep/2,w,g,q)
let xk4=tstep*(v+vk3)
let vk4=tstep*accel(x+xk3,v+vk3,t+tstep,w,g,q)
let vnew=v+(vk1+2*vk2+2*vk3+vk4)/6
let xnew=x+(xk1+2*xk2+2*xk3+xk4)/6
end sub
!
!
def accel(x,v,t,w,g,q)
let accel=-sin(x)-(1/q)*v+g*cos(w*t)
end def